Introduction to Solving Biological Problems Using R - Day 1

Mark Dunning and Aiora Zabala. Original material by Robert Stojnić, Laurent Gatto, Rob Foy John Davey, Dávid Molnár and Ian Roberts

Last modified: 12 Jun 2015

Day 1 Schedule

Introduction to R and its environment
Data Structures
Data Analysis Example
Plotting in R

What’s R?

A statistical programming environment
- based on ‘S’
- suited to high-level data analysis
Open source and cross platform
Extensive graphics capabilities
Diverse range of add-on packages
Active community of developers
Thorough documentation

The R-project page

http://www.r-project.org/

R screenshot

R in the New York Times

http://goo.gl/pww4ZO

R new york time

R in Nature

R in nature

R plotting capabilities

http://spatial.ly/2012/02/great-maps-ggplot2/ R cycle journeys

R plotting capabilities

https://www.facebook.com/notes/facebook-engineering/visualizing-friendships/469716398919 R facebook

Who uses R? Not just academics!

http://www.revolutionanalytics.com/companies-using-r

Various platforms supported

Release 3.2.0 (April 2015)
- Base package and Contributed packages (general purpose extras)
- 6770 available packages as of Fri Jun 12 16:10:52 2015
Download from http://mirrors.ebi.ac.uk/CRAN/
Windows, Mac and Linux versions available
Executed using command line, or a graphical user interface (GUI)
On this course, we use the RStudio GUI (www.rstudio.com)
Everything you need is installed on the training machines
If you are using your own machine, download both R and RStudio

Getting started

R is a program which, once installed on your system, can be launched and is immediately ready to take input directly from the user
There are two ways to launch R:
- From the command line (particularly useful if you’re quite familiar with Linux; in the console at the prompt simply type R)
- As an application called (very good for beginners)

Launching R Using RStudio

To launch RStudio, find the RStudio icon in the menu bar on the left of the screen and click

RStudio screenshot

The Working Directory (wd)

Like many programs R has a concept of a working directory (wd)
It is the place where R will look for files to execute and where it will save files, by default
For this course we need to set the working directory to the location of the course scripts
At the command prompt in the terminal or in RStudio console type:

setwd("R_course/Day_1_scripts")

Alternatively in RStudio use the mouse and browse to the directory location
Session → Set Working Directory → Choose Directory…

Basic concepts in R - command line calculation

The command line can be used as a calculator. Type:

2 + 2

## [1] 4

20/5 - sqrt(25) + 3^2

## [1] 8

sin(pi/2)

## [1] 1

Note: The number in the square brackets is an indicator of the position in the output. In this case the output is a ‘vector’ of length 1 (i.e. a single number). More on vectors coming up…

Basic concepts in R - variables

A variable is a letter or word which takes (or contains) a value. We use the assignment ‘operator’, <-

x <- 10
x

## [1] 10

myNumber <- 25
myNumber

## [1] 25

We can perform arithmetic on variables:

sqrt(myNumber)

## [1] 5

We can add variables together:

x + myNumber

## [1] 35

Basic concepts in R - variables

We can change the value of an existing variable:

x <- 21
x

## [1] 21

We can set one variable to equal the value of another variable:

x <- myNumber
x

## [1] 25

We can modify the contents of a variable:

myNumber <- myNumber + sqrt(16)
myNumber

## [1] 29

Basic concepts in R - functions

Functions in R perform operations on arguments (the inputs(s) to the function). We have already used:

sin(x)

this returns the sine of x. In this case the function has one argument: x. Arguments are always contained in parentheses – curved brackets, () – separated by commas.

Try these:

sum(3,4,5,6)

## [1] 18

max(3,4,5,6)

## [1] 6

min(3,4,5,6)

## [1] 3

Basic concepts in R - functions

Arguments can be named or unnamed, but if they are unnamed they must be ordered (we will see later how to find the right order)

seq(from = 2, to = 20, by = 4)

## [1]  2  6 10 14 18

seq(2, 20, 4)

## [1]  2  6 10 14 18

Basic concepts in R - vectors

The basic data structure in R is a vector – an ordered collection of values. R treats even single values as 1-element vectors. The function c combines its arguments into a vector:

x <- c(3,4,5,6)
x

## [1] 3 4 5 6

The square brackets [] indicate the position within the vector (the index). We can extract individual elements by using the [] notation:

x[1]

## [1] 3

x[4]

## [1] 6

We can even put a vector inside the square brackets: (vector indexing)

y <- c(2,3)
x[y]

## [1] 4 5

Basic concepts in R - vectors

There are a number of shortcuts to create a vector. Instead of:

x <- c(3,4,5,6,7,8,9,10,11,12)

we can write:

x <- 3:12
x

##  [1]  3  4  5  6  7  8  9 10 11 12

or we can use the seq() function, which returns a vector:

x <- seq(2, 20, 4)
x

## [1]  2  6 10 14 18

x <- seq(2, 20, length.out=5)
x

## [1]  2.0  6.5 11.0 15.5 20.0

Basic concepts in R - vectors

or we can use the rep() function:

y <- rep(3, 5)
y

## [1] 3 3 3 3 3

y <- rep(1:3, 5)
y

##  [1] 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3

Basic concepts in R - vectors

We have seen some ways of extracting elements of a vector. We can use these shortcuts to make things easier (or more complex!)

x <- 3:12
x[3:7]

## [1] 5 6 7 8 9

x[seq(2, 6, 2)]

## [1] 4 6 8

x[rep(3, 2)]

## [1] 5 5

Basic concepts in R - vectors

We can add an element to a vector:

y <- c(x, 1)
y

##  [1]  3  4  5  6  7  8  9 10 11 12  1

We can glue vectors together:

z <- c(x, y)
z

##  [1]  3  4  5  6  7  8  9 10 11 12  3  4  5  6  7  8  9 10
## [19] 11 12  1

Basic concepts in R - vectors

We can remove element(s) from a vector:

x <- 3:12
x[-3]

## [1]  3  4  6  7  8  9 10 11 12

x[-(5:7)]

## [1]  3  4  5  6 10 11 12

x[-seq(2, 6, 2)]

## [1]  3  5  7  9 10 11 12

Basic concepts in R - vectors

Finally, we can modify the contents of a vector:

x[6] <- 4
x

##  [1]  3  4  5  6  7  4  9 10 11 12

x[3:5] <- 1
x

##  [1]  3  4  1  1  1  4  9 10 11 12

Remember! Square brackets for indexing [], parentheses for function arguments ().

Basic concepts in R - vector arithmetic

When applying all standard arithmetic operations to vectors, application is element-wise

x <- 1:10
y <- x*2
y

##  [1]  2  4  6  8 10 12 14 16 18 20

z <- x^2
z

##  [1]   1   4   9  16  25  36  49  64  81 100

Basic concepts in R - vector arithmetic

Adding two vectors

y + z

##  [1]   3   8  15  24  35  48  63  80  99 120

If vectors are not the same length, the shorter one will be recycled:

x + 1:2

##  [1]  2  4  4  6  6  8  8 10 10 12

But be careful if the vector lengths aren’t factors of each other:

x + 1:3

## Warning in x + 1:3: longer object length is not a
## multiple of shorter object length

##  [1]  2  4  6  5  7  9  8 10 12 11

Basic concepts in R - Character vectors and naming

All the vectors we have seen so far have contained numbers, but we can also store strings in vectors – this is called a character vector.

gene.names <- c("Pax6", "Beta-actin", "FoxP2", "Hox9")

We can name elements of vectors using the names function, which can be useful to keep track of the meaning of our data:

gene.expression <- c(0, 3.2, 1.2, -2)
gene.expression

## [1]  0.0  3.2  1.2 -2.0

names(gene.expression) <- gene.names
gene.expression

##       Pax6 Beta-actin      FoxP2       Hox9 
##        0.0        3.2        1.2       -2.0

We can also use the names function to get a vector of the names of an object:

names(gene.expression)

## [1] "Pax6"       "Beta-actin" "FoxP2"     
## [4] "Hox9"

Exercise: genes and genomes

Let’s try some vector arithmetic. Here are the genome lengths and number of protein coding genes for several model organisms:

Species	Genome size (Mb)	Protein coding genes
Homo sapiens	3,102	20,774
Mus musculus	2,731	23,139
Drosophila melanogaster	169	13,937
Caenorhabditis elegans	100	20,532
Saccharomyces cerevisiae	12	6,692

Create genome.size and coding.genes vectors to hold the data in each column using the c function. Create a species.name vector and use this vector to name the values in the other two vectors.

Exercise: genes and genomes

Let’s assume a coding gene has an average length of 1.5 kilobases. On average, how many base pairs of each genome is made of coding genes? Create a new vector to record this called coding.bases.
What percentage of each genome is made up of protein coding genes? Use your coding.bases and genome.size vectors to calculate this. (See earlier slides for how to do division in R.)
How many times more bases are used for coding in the human genome compared to the yeast genome? How many times more bases are in the human genome in total compared to the yeast genome? Look up indices of your vectors to find out.

Answers to genome exercise

genome.size  <- c(3102, 2731, 169, 100, 12)
coding.genes <- c(20774, 23139, 13937, 20532, 6692)
species.name <- c("H. sapiens","M. musculus",
                  "D. melanogaster","C. elegans",
                  "S. cerevisiae")
names(genome.size)  <- species.name
names(coding.genes) <- species.name

To calculate the number of coding bases, we need to use the same scale as we used for genome size: 1.5 kilobases is 0.0015 Megabases.

coding.bases <- coding.genes*0.0015
coding.bases

##      H. sapiens     M. musculus D. melanogaster 
##         31.1610         34.7085         20.9055 
##      C. elegans   S. cerevisiae 
##         30.7980         10.0380

Answers to genome exercise

To calculate the percentage of coding bases in each genome:

coding.pc <- coding.bases/genome.size*100
coding.pc

##      H. sapiens     M. musculus D. melanogaster 
##        1.004545        1.270908       12.370118 
##      C. elegans   S. cerevisiae 
##       30.798000       83.650000

To compare human to yeast:

coding.bases[1]/coding.bases[5]

## H. sapiens 
##   3.104304

genome.size[1]/genome.size[5]

## H. sapiens 
##      258.5

Note that if a new vector is created using a named vector, the names are usually carried across to the new vector. Sometimes this is what we want (as for coding.pc) but sometimes it is not (when we are comparing human to yeast). We can remove names by setting them to the special NULL value:

names(coding.pc) <- NULL
coding.pc

## [1]  1.004545  1.270908 12.370118 30.798000
## [5] 83.650000

Writing scripts with RStudio

Typing lots of commands directly to R can be tedious. A better way is to write the commands to a file and then load it into R.

Click on File → New in Rstudio
Type in some R code, e.g.:

x <- 2 + 2
print(x)

Click on Run to execute the current line, and Source to execute the whole script

Sourcing can also be performed manually with source("myScript.R")

Getting help

This is possibly the most important slide in the whole course
To get help on any R function, type ? followed by the function name. For example:

?seq

This retrieves the syntax and arguments for the function. You can see the default order of arguments here. The help page also tells you which package it belongs to.
There will typically be example usage, which you can test using the example function:

example(seq)

If you can’t remember the exact name type ?? followed by your guess. R will return a list of possibilities:

??plot

The Packages tab in the lower-right panel of RStudio will help you to locate the help pages for a particular package and its functions
- Often there will be a user-guide / ‘vignette’ too

Interacting with the R console

R console symbols:
- ; end of line (Enables multiple commands to be placed on one line of text)
- # comment (indicates text is a comment and not executed)
- + command line wrap (R is waiting for you to complete an expression)
Ctrl-c or escape to clear input line and try again
Ctrl-l to clear window
Use the TAB key for command auto completion
Use up and down arrows to scroll through the command history

R packages

R comes ready loaded with various libraries of functions called packages. e.g. the function sum() is in the base package and sd(), which calculates the standard deviation of a vector, is in the stats package
There are 1000s of additional packages provided by third parties, and the packages can be found in numerous server locations on the web called repositories
The two repositories you will come across the most are
- The Comprehensive R Archive Network (CRAN)
- Bioconductor
CRAN packages can be installed using install.packages
Set the Bioconductor package download tool by typing:

source("http://bioconductor.org/biocLite.R")

Bioconductor packages are then installed with the biocLite() function:

biocLite("PackageName")

Exercise: Install packages ggplot2 and DESeq

Packages on CRAN:
- Use CRAN search to find functionality you need: http://cran.r-project.org/search.html
- Or look for packages by theme: http://cran.r-project.org/web/views/
Packages in Bioconductor – specialised in genomics: http://www.bioconductor.org/packages/release/bioc/
Other repositories:
- http://r-forge.r-project.org/
- github.com can also host R packages
Bottomline: always first look if there is already an R package that does what you want before trying to implement it yourself

Exercise: Install packages ggplot2 and DESeq

ggplot2 is a commonly used graphics package
- in RStudio, go to Tools → Install Packages… and type the package name
- or use install.packages() function to install it:

install.packages("ggplot2")

DESeq is a Bioconductor package (http://www.bioconductor.org) for the analysis of RNA-seq data

source("http://www.bioconductor.org/biocLite.R")
biocLite("DESeq")

R needs to be told to use the new functions from the installed packages. Use library(...) function to load the newly installed features:

library(ggplot2) # loads ggplot functions
library(DESeq)   # loads DESeq functions
library()        # Lists all the packages you've got installed

R is designed to handle experimental data

Although the basic unit of R is a vector, we usually handle data in data frames.
A data frame is a set of observations of a set of variables – in other words, the outcome of an experiment.
For example, we might want to analyse information about a set of patients. To start with, let’s say we have ten patients and for each one we know their name, sex, age, weight and whether they give consent for their data to be made public.

The patients data frame

We are going to create a data frame called ‘patients’, which will have ten rows (observations) and seven columns (variables). The columns must all be equal lengths.
We will explore how to construct these data from scratch.

	First_Name	Second_Name	Full_Name	Sex	Age	Weight	Consent
1	Adam	Jones	Adam Jones	Male	50	70.8	TRUE
2	Eve	Parker	Eve Parker	Female	21	67.9	TRUE
3	John	Evans	John Evans	Male	35	75.3	FALSE
4	Mary	Davis	Mary Davis	Female	45	61.9	TRUE
5	Peter	Baker	Peter Baker	Male	28	72.4	FALSE
6	Paul	Daniels	Paul Daniels	Male	31	69.9	FALSE
7	Joanna	Edwards	Joanna Edwards	Female	42	63.5	FALSE
8	Matthew	Smith	Matthew Smith	Male	33	71.5	TRUE
9	David	Roberts	David Roberts	Male	57	73.2	FALSE
10	Sally	Wilson	Sally Wilson	Female	62	64.8	TRUE

Character, numeric and logical data types

Each column is a vector, like previous vectors we have seen, for example:

age    <- c(50, 21, 35, 45, 28, 31, 42, 33, 57, 62)
weight <- c(70.8, 67.9, 75.3, 61.9, 72.4, 69.9, 63.5, 
71.5, 73.2, 64.8)

We can define the names using character vectors:

firstName <- c("Adam", "Eve", "John", "Mary", "Peter", 
"Paul", "Joanna", "Matthew", "David", "Sally")
secondName <- c("Jones", "Parker", "Evans", "Davis", 
"Baker","Daniels", "Edwards", "Smith", "Roberts", "Wilson")

We also have a new type of vector, the logical vector, which only contains the values TRUE and FALSE:

consent <- c(TRUE, TRUE, FALSE, TRUE, FALSE, FALSE,
FALSE, TRUE, FALSE, TRUE)

Character, numeric and logical data types

Vectors can only contain one type of data; we cannot mix numbers, characters and logical values in the same vector. If we try this, R will convert everything to characters:

c(20, "a string", TRUE)

## [1] "20"       "a string" "TRUE"

We can see the type of a particular vector using the class function:

 class(firstName)

## [1] "character"

 class(age)

## [1] "numeric"

 class(weight)

## [1] "numeric"

 class(consent)

## [1] "logical"

Factors

Character vectors are fine for some variables, like names
But sometimes we have categorical data and we want R to recognize this.
A factor is R’s data structure for categorical data.

sex <- c("Male", "Female", "Male", "Female", "Male", "Male",
"Female", "Male", "Male", "Female")

sex

##  [1] "Male"   "Female" "Male"   "Female" "Male"   "Male"  
##  [7] "Female" "Male"   "Male"   "Female"

factor(sex)

##  [1] Male   Female Male   Female Male   Male   Female Male  
##  [9] Male   Female
## Levels: Female Male

R has converted the strings of the sex character vector into two levels, which are the categories in the data
Note the values of this factor are not character strings, but levels
We can use this factor to compare data for males and females

Creating a data frame (first attempt)

We can construct a data frame from other objects (N.B. The paste function joins character vectors together)

patients <- data.frame(firstName, secondName, 
                       paste(firstName, secondName),  
                       sex, age, weight, consent)
patients

##    firstName secondName paste.firstName..secondName.    sex
## 1       Adam      Jones                   Adam Jones   Male
## 2        Eve     Parker                   Eve Parker Female
## 3       John      Evans                   John Evans   Male
## 4       Mary      Davis                   Mary Davis Female
## 5      Peter      Baker                  Peter Baker   Male
## 6       Paul    Daniels                 Paul Daniels   Male
## 7     Joanna    Edwards               Joanna Edwards Female
## 8    Matthew      Smith                Matthew Smith   Male
## 9      David    Roberts                David Roberts   Male
## 10     Sally     Wilson                 Sally Wilson Female
##    age weight consent
## 1   50   70.8    TRUE
## 2   21   67.9    TRUE
## 3   35   75.3   FALSE
## 4   45   61.9    TRUE
## 5   28   72.4   FALSE
## 6   31   69.9   FALSE
## 7   42   63.5   FALSE
## 8   33   71.5    TRUE
## 9   57   73.2   FALSE
## 10  62   64.8    TRUE

Naming data frame variables

We can access particular variables using the ‘$’ operator:

patients$age

##  [1] 50 21 35 45 28 31 42 33 57 62

R has inferred the names of our data frame variables from the names of the vectors or the commands (e.g. the paste command)
We can name the variables after we have created a data frame using the names function, and we can use the same function to see the names:

names(patients) <- c("First_Name", "Second_Name", 
"Full_Name", "Sex", "Age", "Weight", "Consent")

names(patients)

## [1] "First_Name"  "Second_Name"
## [3] "Full_Name"   "Sex"        
## [5] "Age"         "Weight"     
## [7] "Consent"

Naming data frame variables

Or we can name the variables when we define the data frame

patients <- data.frame(First_Name = firstName, 
                       Second_Name = secondName, 
                       Full_Name = paste(firstName,secondName), 
                       Sex = sex,
                       Age = age,
                       Weight = weight, 
                       Consent = consent)

names(patients)

## [1] "First_Name"  "Second_Name"
## [3] "Full_Name"   "Sex"        
## [5] "Age"         "Weight"     
## [7] "Consent"

Factors in data frames

When creating a data frame, R assumes all character vectors should be categorical variables and converts them to factors. This is not always what we want:
- e.g. we are unlikely to be interested in the hypothesis that people called Adam are taller, so it seems a bit silly to represent this as a factor

patients$First_Name

##  [1] Adam    Eve     John    Mary   
##  [5] Peter   Paul    Joanna  Matthew
##  [9] David   Sally  
## 10 Levels: Adam David Eve ... Sally

Factors in data frames

We can avoid this by asking R not to treat strings as factors, and then explicitly stating when we want a factor by using factor:

patients <- data.frame(First_Name = firstName, 
                       Second_Name = secondName, 
                       Full_Name = paste(firstName, secondName), 
                       Sex = factor(sex),
                       Age = age,
                       Weight = weight,
                       Consent = consent,
                       stringsAsFactors = FALSE)

 patients$Sex

##  [1] Male   Female Male   Female Male  
##  [6] Male   Female Male   Male   Female
## Levels: Female Male

 patients$First_Name

##  [1] "Adam"    "Eve"     "John"   
##  [4] "Mary"    "Peter"   "Paul"   
##  [7] "Joanna"  "Matthew" "David"  
## [10] "Sally"

Matrices

Data frames are R’s speciality, but R also handles matrices:
- all columns are assumed to contain the same data type. e.g. numerical
- some calculations are more efficient to do on matrices

e <- matrix(1:10, nrow=5, ncol=2)
e

##      [,1] [,2]
## [1,]    1    6
## [2,]    2    7
## [3,]    3    8
## [4,]    4    9
## [5,]    5   10

f <- matrix(1:10, nrow=2, ncol=5)
f

##      [,1] [,2] [,3] [,4] [,5]
## [1,]    1    3    5    7    9
## [2,]    2    4    6    8   10

f %*% e

##      [,1] [,2]
## [1,]   95  220
## [2,]  110  260

The %*% operator is the matrix multiplication operator, not the standard multiplication operator.

Lists

We have seen that vectors can only hold data of one type. How can we store data of multiple types? Or vectors of different lengths in one object?
We can use lists. A list can contain objects of any type

one.to.ten <- 1:10
uniform.mat <- matrix(runif(100), ncol=10, nrow=10)
year.to.october <- data.frame(one.to.ten, month.name[1:10])
myList <- list(ls.obj.1 = one.to.ten, 
               ls.obj.2 = uniform.mat, 
               ls.obj.3 = year.to.october)
names(myList)

## [1] "ls.obj.1" "ls.obj.2" "ls.obj.3"

Lists

We can use the dollar syntax to access list items (in fact, a data frame is a special type of list):

myList$ls.obj.1

##  [1]  1  2  3  4  5  6  7  8  9 10

We can also use myList[[1]] to get the first item in the list.

 myList[[1]]

##  [1]  1  2  3  4  5  6  7  8  9 10

Indexing data frames and matrices

You can index multidimensional data structures like matrices and data frames using commas. If you don’t provide an index for either rows or columns, all of the rows or columns will be returned.

object[rows, colums]

e[1,2]

## [1] 6

e[1,]

## [1] 1 6

patients[1,2]

## [1] "Jones"

patients[1,]

##   First_Name Second_Name  Full_Name
## 1       Adam       Jones Adam Jones
##    Sex Age Weight Consent
## 1 Male  50   70.8    TRUE

Advanced indexing

As values in R are really vectors, so indices are actually vectors, and can be numeric or logical:

 s <- letters[1:5]
 s[c(1,3)]

## [1] "a" "c"

 s[c(TRUE, FALSE, TRUE, FALSE, FALSE)]

## [1] "a" "c"

a <- 1:5
a < 3

## [1]  TRUE  TRUE FALSE FALSE FALSE

s[a < 3]

## [1] "a" "b"

s[a > 1 & a <3]

## [1] "b"

s[a == 2]

## [1] "b"

Operators

arithmetic +, -, *, /, ^
comparison <, >, <=, >=, ==, !=
logical !, &, |, xor
- The operators for ‘comparison’ and ‘logical’ always return logical values! i.e. (TRUE, FALSE)

Exercise

Create a data frame called my.patients using the instructions in the slides. Change the data if you like.
Check you have created the data frame correctly by loading the original version from this file in the Day_1_scripts folder using source

source("1.2_patients.R")

Remake your data frame with three new variables: country, continent, and height
- Make up the data
- Make country a character vector but continent a factor
Try the summary function on your data frame. What does it do? How does it treat vectors (numeric, character, logical) and factors? (What does it do for matrices?)
Use logical indexing to select the following patients:
- Patients under 40
- Patients who give consent to share their data
- Men who weigh as much or more than the average European male (70.8 kg)

Logical indexing answers

Patients under 40:

patients[patients$Age < 40,]

Patients who give consent to share their data:

patients[patients$Consent == TRUE,]

Men who weigh as much or more than the average European male (70.8 kg):

patients[patients$Sex == "Male" & patients$Weight >= 70.8,]

Basic R ‘Built-in’ functions for working with data frames

We have seen before how we can get the names of our variables, but for data frames and matrices we can also get these names with colnames, and the names of the rows with rownames

names(patients)

## [1] "First_Name"  "Second_Name"
## [3] "Full_Name"   "Sex"        
## [5] "Age"         "Weight"     
## [7] "Consent"

colnames(patients)

## [1] "First_Name"  "Second_Name"
## [3] "Full_Name"   "Sex"        
## [5] "Age"         "Weight"     
## [7] "Consent"

rownames(patients)

##  [1] "1"  "2"  "3"  "4"  "5"  "6"  "7" 
##  [8] "8"  "9"  "10"

Basic R ‘Built-in’ functions for working with data frames

We can get the numbers of rows or columns with nrow and ncol

nrow(patients)

## [1] 10

ncol(patients)

## [1] 7

dim(patients)

## [1] 10  7

Basic R ‘Built-in’ functions for working with data frames

We can also find the length of a vector or a list with length, although this may give confusing results for some structures, like data frames:

length(c(1,2,3,4,5))

## [1] 5

length(patients)

## [1] 7

length(patients$Age)

## [1] 10

Remember, a data frame is a list of variables, so its length is the number of variables. The length of one of the variable vectors (like Age) is the number of observations

Basic R ‘Built-in’ functions for working with data frames

We can add rows or columns to a data frame using rbind and cbind

newpatient <- c("Kate", "Lawson", "Kate Lawson", "Female", "35", "62.5","TRUE")
rbind(patients,newpatient)

##    First_Name Second_Name      Full_Name    Sex Age Weight Consent
## 1        Adam       Jones     Adam Jones   Male  50   70.8    TRUE
## 2         Eve      Parker     Eve Parker Female  21   67.9    TRUE
## 3        John       Evans     John Evans   Male  35   75.3   FALSE
## 4        Mary       Davis     Mary Davis Female  45   61.9    TRUE
## 5       Peter       Baker    Peter Baker   Male  28   72.4   FALSE
## 6        Paul     Daniels   Paul Daniels   Male  31   69.9   FALSE
## 7      Joanna     Edwards Joanna Edwards Female  42   63.5   FALSE
## 8     Matthew       Smith  Matthew Smith   Male  33   71.5    TRUE
## 9       David     Roberts  David Roberts   Male  57   73.2   FALSE
## 10      Sally      Wilson   Sally Wilson Female  62   64.8    TRUE
## 11       Kate      Lawson    Kate Lawson Female  35   62.5    TRUE

Note we are not altering the patients data frame, just displaying to the screen
- we would have to use an assignment <- operator to create a new object

Basic R ‘Built-in’ functions for working with data frames

cbind(patients,10:1)

##    First_Name Second_Name      Full_Name    Sex Age Weight Consent 10:1
## 1        Adam       Jones     Adam Jones   Male  50   70.8    TRUE   10
## 2         Eve      Parker     Eve Parker Female  21   67.9    TRUE    9
## 3        John       Evans     John Evans   Male  35   75.3   FALSE    8
## 4        Mary       Davis     Mary Davis Female  45   61.9    TRUE    7
## 5       Peter       Baker    Peter Baker   Male  28   72.4   FALSE    6
## 6        Paul     Daniels   Paul Daniels   Male  31   69.9   FALSE    5
## 7      Joanna     Edwards Joanna Edwards Female  42   63.5   FALSE    4
## 8     Matthew       Smith  Matthew Smith   Male  33   71.5    TRUE    3
## 9       David     Roberts  David Roberts   Male  57   73.2   FALSE    2
## 10      Sally      Wilson   Sally Wilson Female  62   64.8    TRUE    1

Note we are not altering the patients data frame, just displaying to the screen
- we would have to use an assignment <- operator to create a new object

Basic R ‘Built-in’ functions for working with data frames

We can also remove rows and columns

patients[-1, ] # remove first row

##    First_Name Second_Name      Full_Name    Sex Age Weight Consent
## 2         Eve      Parker     Eve Parker Female  21   67.9    TRUE
## 3        John       Evans     John Evans   Male  35   75.3   FALSE
## 4        Mary       Davis     Mary Davis Female  45   61.9    TRUE
## 5       Peter       Baker    Peter Baker   Male  28   72.4   FALSE
## 6        Paul     Daniels   Paul Daniels   Male  31   69.9   FALSE
## 7      Joanna     Edwards Joanna Edwards Female  42   63.5   FALSE
## 8     Matthew       Smith  Matthew Smith   Male  33   71.5    TRUE
## 9       David     Roberts  David Roberts   Male  57   73.2   FALSE
## 10      Sally      Wilson   Sally Wilson Female  62   64.8    TRUE

Basic R ‘Built-in’ functions for working with data frames

We can also remove rows and columns

patients[ ,-1] # remove first column

##    Second_Name      Full_Name    Sex Age Weight Consent
## 1        Jones     Adam Jones   Male  50   70.8    TRUE
## 2       Parker     Eve Parker Female  21   67.9    TRUE
## 3        Evans     John Evans   Male  35   75.3   FALSE
## 4        Davis     Mary Davis Female  45   61.9    TRUE
## 5        Baker    Peter Baker   Male  28   72.4   FALSE
## 6      Daniels   Paul Daniels   Male  31   69.9   FALSE
## 7      Edwards Joanna Edwards Female  42   63.5   FALSE
## 8        Smith  Matthew Smith   Male  33   71.5    TRUE
## 9      Roberts  David Roberts   Male  57   73.2   FALSE
## 10      Wilson   Sally Wilson Female  62   64.8    TRUE

Basic R ‘Built-in’ functions for working with data frames

Sorting a vector with sort

sort(patients$Second_Name)

##  [1] "Baker"   "Daniels" "Davis"   "Edwards" "Evans"   "Jones"   "Parker" 
##  [8] "Roberts" "Smith"   "Wilson"

Sorting a data frame by one variable with order

order(patients$Second_Name)

##  [1]  5  6  4  7  3  1  2  9  8 10

patients[order(patients$Second_Name),]

##    First_Name Second_Name      Full_Name    Sex Age Weight Consent
## 5       Peter       Baker    Peter Baker   Male  28   72.4   FALSE
## 6        Paul     Daniels   Paul Daniels   Male  31   69.9   FALSE
## 4        Mary       Davis     Mary Davis Female  45   61.9    TRUE
## 7      Joanna     Edwards Joanna Edwards Female  42   63.5   FALSE
## 3        John       Evans     John Evans   Male  35   75.3   FALSE
## 1        Adam       Jones     Adam Jones   Male  50   70.8    TRUE
## 2         Eve      Parker     Eve Parker Female  21   67.9    TRUE
## 9       David     Roberts  David Roberts   Male  57   73.2   FALSE
## 8     Matthew       Smith  Matthew Smith   Male  33   71.5    TRUE
## 10      Sally      Wilson   Sally Wilson Female  62   64.8    TRUE

Exercise

Output the patients data frame, with the patients sorted in order of age, oldest first. (You may need the rev function, or look at the options for the order function)
Add a height column to the data frame using the cbind function and assign the result to a new data frame
- enter your own height values, or use these 1.8, 1.6, 1.9, 1.5, 1.7, 1.5,1.75,1.8,1.6,1.5
Calculate the body mass index for each patient. Order the patients according to their bmi
- $bmi = \frac{weight}{height^2}$
A person is classified as overweight if they have a bmi > 25. Which patients are overweight using this criteria, and have consented to having their data shared?

Answers

Make sure we have the correct patients data frame

source("1.2_patients.R")

patients[order(patients$Age,decreasing = TRUE),]

##    First_Name Second_Name      Full_Name    Sex Age Weight Consent
## 10      Sally      Wilson   Sally Wilson Female  62   64.8    TRUE
## 9       David     Roberts  David Roberts   Male  57   73.2   FALSE
## 1        Adam       Jones     Adam Jones   Male  50   70.8    TRUE
## 4        Mary       Davis     Mary Davis Female  45   61.9    TRUE
## 7      Joanna     Edwards Joanna Edwards Female  42   63.5   FALSE
## 3        John       Evans     John Evans   Male  35   75.3   FALSE
## 8     Matthew       Smith  Matthew Smith   Male  33   71.5    TRUE
## 6        Paul     Daniels   Paul Daniels   Male  31   69.9   FALSE
## 5       Peter       Baker    Peter Baker   Male  28   72.4   FALSE
## 2         Eve      Parker     Eve Parker Female  21   67.9    TRUE

patients2 <- cbind(patients, Height = c(1.8, 1.6, 1.9, 1.5, 1.7, 1.5,1.75,1.8,1.6,1.5))
patients2

##    First_Name Second_Name      Full_Name    Sex Age Weight Consent Height
## 1        Adam       Jones     Adam Jones   Male  50   70.8    TRUE   1.80
## 2         Eve      Parker     Eve Parker Female  21   67.9    TRUE   1.60
## 3        John       Evans     John Evans   Male  35   75.3   FALSE   1.90
## 4        Mary       Davis     Mary Davis Female  45   61.9    TRUE   1.50
## 5       Peter       Baker    Peter Baker   Male  28   72.4   FALSE   1.70
## 6        Paul     Daniels   Paul Daniels   Male  31   69.9   FALSE   1.50
## 7      Joanna     Edwards Joanna Edwards Female  42   63.5   FALSE   1.75
## 8     Matthew       Smith  Matthew Smith   Male  33   71.5    TRUE   1.80
## 9       David     Roberts  David Roberts   Male  57   73.2   FALSE   1.60
## 10      Sally      Wilson   Sally Wilson Female  62   64.8    TRUE   1.50

Answers

bmi <- patients2$Weight / (patients2$Height)^2
bmi

##  [1] 21.85185 26.52344 20.85873 27.51111 25.05190 31.06667 20.73469 22.06790
##  [9] 28.59375 28.80000

patients2[bmi > 25 & patients2$Consent,]

##    First_Name Second_Name    Full_Name    Sex Age Weight Consent Height
## 2         Eve      Parker   Eve Parker Female  21   67.9    TRUE    1.6
## 4        Mary       Davis   Mary Davis Female  45   61.9    TRUE    1.5
## 10      Sally      Wilson Sally Wilson Female  62   64.8    TRUE    1.5

3 steps to Basic Data Analysis

In this short section, we show how the data manipulation steps we have just seen can be used as part of an analysis pipeline

Reading in data
- read.table()
- read.csv(), read.delim()
Analysis
- Manipulating & reshaping the data
- Any maths you like
- Plotting the outcome
Writing out results
- write.table()
- write.csv()

A simple walkthrough

50 neuroblastoma patients were tested for NMYC gene copy number by interphase nuclei FISH
- Amplification of NMYC correlates with worse prognosis
- We have count data
- Numbers of cells per patient assayed
  - For each we have NMYC copy number relative to base ploidy
We need to determine which patients have amplifications
- (i.e > 33% of cells show NMYC amplification)

1. Read in the data

The data is a tab-delimited file. Each row is a record, each column is a field. Columns are separated by tabs in the text
We need to read in the results and assign it to an object (rawdata)

rawData <- read.delim("1.3_NBcountData.txt")
rawData[1:10,]

If the data has been comma-separated then, sep=",":

read.csv("1.3_NBcountData.csv")

For full list of arguments

?read.table

1. Read in the data

rawData[1:10,]  # View the first 10 rows to ensure import is OK

##    Patient Nuclei NB_Amp NB_Nor NB_Del
## 1        1     65      0     63      2
## 2        2     51      3     43      5
## 3        3     37      2     35      0
## 4        4     37      2     35      0
## 5        5     45      2     42      1
## 6        6     46      4     41      1
## 7        7     65      1     64      0
## 8        8     59      1     54      4
## 9        9     49      0     48      1
## 10      10     46      0     45      1

Word of caution

tolstoy

hadley

Like families, tidy datasets are all alike but every messy dataset is messy in its own way - (Hadley Wickham)

http://vimeo.com/33727555

http://vita.had.co.nz/papers/tidy-data.pdf

Handling missing values

The data frame contains some NA values, which means the values are missing – a common occurrence in real data collection
NA is a special value that can be present in objects of any type (logical, character, numeric etc)
NA is not the same as NULL. NULL is an empty R object. NA is one missing value within an R object (like a data frame or a vector)
Often R functions will handle NAs gracefully, but sometimes we have to tell the functions what to do with them. R has some built-in functions for dealing with NAs, and functions often have their own arguments (like na.rm) for handling them

x <- c(1,NA,3)
mean(x)

## [1] NA

mean(x, na.rm=TRUE)

## [1] 2

mean(na.omit(x))

## [1] 2

 is.na(x)

## [1] FALSE  TRUE FALSE

2. Analysis (reshaping data and maths)

Our analysis involves identifying patients with > 33% NB amplification

# create an index of results:
prop <- rawData$NB_Amp / rawData$Nuclei 
# Get sample names of amplified patients:
amp <- which(prop > 0.33)

We can plot a simple chart of the % NB amplification
- note that two samples are amplified
- plotting will be covered in detail shortly

plot(prop, ylim=c(0,1))
abline(h=0.33) # Add a horizonal line

3. Outputting the results

We write out a data frame of results (patients > 33% NB amplification) as a ‘comma separated values’ text file (CSV):

write.csv(rawData[amp,], file="selectedSamples.csv")

The output file is directly-readable by Excel
It’s often helpful to double check where the data has been saved.
- Use the get working directory function:

getwd() # print working directory
list.files() # list files in working directory

Data analysis exercise: Which samples are near normal?

Patients are near normal if: (NB_Amp / Nuclei < 0.33 & NB_Del == 0)
Modify the condition in our previous code to find these patients
Write out a results file of the samples that match these criteria, and open it in a spreadsheet program

Solution to NB normality test

norm <- which(prop < 0.33 & rawData$NB_Del == 0)
norm

##  [1]  3  4  7 15 20 24 36 37 42 47

write.csv(rawData[norm,], "My_NB_output.csv")

Plot basics

As we have heard, R has extensive graphical capabilities
but we need to start simple
we will describe base graphics in R; the plots available with any standard R installation
- other alternatives are available but require more R knowledge; e.g. lattice, ggplot2
plot is the most basic plotting function in R
- it can try and guess the best plot based on the data

Making a Scatter Plot

If given a single vector as an argument, plot will make a scatter plot with the values of the vector on the y axis, and indices in the x axis
- e.g. it puts a point at
  - x=1,y=70.8
  - x=2,y=67.9 etc…

patients$Weight

##  [1] 70.8 67.9 75.3 61.9 72.4 69.9 63.5 71.5 73.2 64.8

plot(patients$Weight)

Making a Scatter Plot

R tries to guess the most appropriate way to visualise the data

Axis limits, labels, titles are inferred from the data
- we can fix these later on

Making a Scatter Plot of two variables

We can give two arguments to plot. It will put the values from the first argument in the x axis, and values from the second argument on the y axis
- for when we want to visualise the relationship between two variables

patients$Age

##  [1] 50 21 35 45 28 31 42 33 57 62

plot(patients$Weight,patients$Age)

Making a barplot

barplot(patients$Age)

Making a barplot

It is more usual to display count data in a barplot
- e.g. the counts of the a particular categorical variable

barplot(summary(patients$Sex))

Plotting a distribution: Histogram

A histogram is a popular way of visualising a distribution

hist(patients$Weight)

Plotting a distribution: Histogram

A histogram is a popular way of visualising a distribution

hist(patients$Age)

Plotting a distribution: Boxplot

The boxplot is commonly-used in statistics to visualise a distribution
- The black solid line is the median
- The top and bottom of the box is the 75th and 25th percentiles
  - Hence, the distance between these is a reflection of the spread of the data; the Inter-Quartile Range (IQR)
- Whiskers are drawn at 1.5 x IQR and -1.5 x IQR
- Any points are outliers

boxplot(patients$Weight)

Plotting a distribution: Boxplot

Sometimes we want to compare distributions between different groupings in our data
For this we need to use the ‘formula’ syntax
- For now, y ~ x means put variable y on the y axis and x on the x axis

boxplot(patients$Weight~patients$Sex)

Exercise: Data exploration

In the Day_1_scripts folder you will find the file ozone.csv
- Data describing weather conditions in New York City in 1973, obtained from the supplementary data to Biostatistics: A Methodology for the Health Sciences
- Full description here http://faculty.washington.edu/heagerty/Books/Biostatistics/DATA/ozonedoc.txt
Load these data into R
What data types are present? Try to think of ways to create the following plots from the data
- scatter plot two variables
- a histogram
- boxplot of a continuous variable against a categorical variable

Suggestions

weather <- read.csv("ozone.csv")
head(weather)

##   Ozone Solar.R Wind Temp Month Day
## 1    41     190  7.4   67     5   1
## 2    36     118  8.0   72     5   2
## 3    12     149 12.6   74     5   3
## 4    18     313 11.5   62     5   4
## 5    NA      NA 14.3   56     5   5
## 6    28      NA 14.9   66     5   6

Suggestions

plot(weather$Ozone,weather$Solar.R)

Suggestions

hist(weather$Temp)

Suggestions

boxplot(weather$Ozone)

Suggestions

boxplot(weather$Ozone~weather$Month)

Simple customisations

plot comes with a whole array of arguments that can be set when we call the function
- see ?plot
we can choose to draw lines on the plot rather than points
- the rest of the plot is the same as before

plot(patients$Weight,type="l")

Simple customisations

We can also have both lines and points

plot(patients$Weight,type="b")

Simple customisations

Add an informative title to the plot using the main argument

plot(patients$Weight,patients$Age,main="Relationship between Weight and Age")

Simple customisations

Adding the x-axis label

plot(patients$Weight,patients$Age,xlab="Weight")

Simple customisations

Adding the y-axis label

plot(patients$Weight,patients$Age,ylab="Age")

Simple customisations

We can specifiy multiple arguments at once
- here ylim and xlim are used to specify axis limits

plot(patients$Weight,patients$Age,
     ylab="Age",
     xlab="Weight",
     main="Relationship between Weight and Age",
     ylim=c(10,70),
     xlim=c(60,80))

Defining a colour

R can recognise various strings "red", "orange","green","blue","yellow"….
Or more exotic ones lightgrey, springgreen4, burlywood3, grey12, tomato, red, palegreen2, hotpink2….. See colours().
See http://www.stat.columbia.edu/~tzheng/files/Rcolor.pdf
Can also use Red Green Blue , hexadecimal, values

Use of colours

Changing the col argument to plot changes the colour that the points are plotted in

plot(patients$Weight,patients$Age,col="red")

Plotting characters

R can use a variety of plotting characters
Each of which has a numeric code

plot(patients$Weight,patients$Age, pch=16)

Plotting characters

Or you can specify a character

plot(patients$Weight,patients$Age, pch="X")

Size of points

Character expansion

plot(patients$Weight,patients$Age, cex=2)

Size of points

Character expansion

plot(patients$Weight,patients$Age, cex=0.2)

Colours and characters as vectors

Previously we have used a vector of length 1 as our value of colour and character
we can use a vector of any length
- the values will get recycled (re-used) so that each point gets assigned a value

plot(patients$Weight,patients$Age, 
     pch = 1:10,
     col=c("red","orange","green","blue"),
     cex = 1:5)

Other plots use the same arguments

Other plotting functions use the same arguments as plot
- technical explanation: the arguments are ‘inherited’

boxplot(patients$Weight~patients$Sex,
        xlab="Sex",
        ylab="Weight",
        main="Relationship between Weight and Gender",
        col=c("blue","yellow"))

Exercise

Take one of the plots you produced in the previous section and experiment with adding axes and colours

Introduction to Solving Biological Problems Using R - Day 1

Day 1 Schedule

1. Introduction to R and its environment

What’s R?

The R-project page

R in the New York Times

R in Nature

R plotting capabilities

R plotting capabilities

Who uses R? Not just academics!

Various platforms supported

Getting started

Launching R Using RStudio

The Working Directory (wd)

Basic concepts in R - command line calculation

Basic concepts in R - variables

Basic concepts in R - variables

Basic concepts in R - functions

Basic concepts in R - functions

Basic concepts in R - vectors

Basic concepts in R - vectors

Basic concepts in R - vectors

Basic concepts in R - vectors

Basic concepts in R - vectors

Basic concepts in R - vectors

Basic concepts in R - vectors

Basic concepts in R - vector arithmetic

Basic concepts in R - vector arithmetic

Basic concepts in R - Character vectors and naming

Exercise: genes and genomes

Exercise: genes and genomes

Answers to genome exercise

Answers to genome exercise

Writing scripts with RStudio

Getting help

Interacting with the R console

R packages

Exercise: Install packages ggplot2 and DESeq

Exercise: Install packages ggplot2 and DESeq

2. Data structures

R is designed to handle experimental data

The patients data frame

Character, numeric and logical data types

Character, numeric and logical data types

Factors

Creating a data frame (first attempt)

Naming data frame variables

Naming data frame variables

Factors in data frames

Factors in data frames

Matrices

Lists

Lists

Indexing data frames and matrices

Advanced indexing

Operators

Exercise

Logical indexing answers

Basic R ‘Built-in’ functions for working with data frames

Basic R ‘Built-in’ functions for working with data frames

Basic R ‘Built-in’ functions for working with data frames

Basic R ‘Built-in’ functions for working with data frames

Basic R ‘Built-in’ functions for working with data frames

Basic R ‘Built-in’ functions for working with data frames

Basic R ‘Built-in’ functions for working with data frames

Basic R ‘Built-in’ functions for working with data frames

Exercise

Answers

Answers

3. R for data analysis

3 steps to Basic Data Analysis

A simple walkthrough

1. Read in the data

1. Read in the data

Word of caution

Handling missing values

2. Analysis (reshaping data and maths)

3. Outputting the results

Data analysis exercise: Which samples are near normal?

Solution to NB normality test